function [ DrN ] = getDrN( r, rn, N, qn )
%GETDRN Summary of this function goes here
%   Detailed explanation goes here

s = 1 / rn;
A = N;
ds = -r / rn^3;
dA2 = [0, dot(qn, r)*r(2)/rn^3, 0, -dot(qn, r)*r(1)/rn^3; -dot(qn, r)*r(2)/rn^3, 0, dot(qn, r)*r(1)/rn^3, 0];
dA = reshape(dA2, 2, 2, 2);

DrN = -1 * tensorProductGradient(s, ds, A, dA);

end

